Chromosome compaction - simple three-state model

<<<<<<< HEAD

library(ChromCom)
library(ggplot2)
library(gridExtra)
library(reshape2)

=======

library(mylib)
library(ChromCom)
library(ggplot2)
library(gridExtra)
library(reshape2)
library(parallel)
binDir <- "../RData"

>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba

11/07/2017

Writing introduction. First script to do simulation.

12/07/2017

Let’s try if it works. An example for \(t_1 = -60\) min, \(\Delta t = 0\), \(k_1 = 0.04\ {\rm min}^{-1}\) and \(k_2 = 0.04\ {\rm min}^{-1}\).

<<<<<<< HEAD

pars <- list(
  t1 = -60,
  dt = 0,
  r1 = 0.04,
  r2 = 0.04
=======

pars <- c3pars(
  t1 = -60,
  dt2 = 0,
  k1 = 0.04,
  k2 = 0.04
>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba
)
chr <- ChromCom3(pars)
chr <- generateCells(chr, nsim=1000)
plotTimelines(chr)


<<<<<<< HEAD

=======

>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba



Just to check if output is as expected. With \(k_2=0\) we should have pure exponential decay. The yellow curve is \(e^{-k_1 t}\).


<<<<<<< HEAD

pars <- list(
  t1 = -60,
  dt = 0,
  r1 = 0.04,
  r2 = 0
=======

pars <- c3pars(
  t1 = -60,
  dt2 = 0,
  k1 = 0.04,
  k2 = 0
>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba
)
tchr <- ChromCom3(pars)
tchr <- generateCells(tchr, nsim=1000)
x <- seq(from=pars$t1, to=100, by=1)
<<<<<<< HEAD
y <- exp(-pars$r1 * (x - pars$t1))
=======
y <- exp(-pars$k1 * (x - pars$t1))
>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba
g <- plotTimelines(tchr)
g + geom_line(data=data.frame(x=x, y=y), aes(x,y), colour="yellow")


<<<<<<< HEAD

=======

>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba



Close, but not perfect. Here is an example with much smaller time step (0.1) and larger number of cells (10,000).


<<<<<<< HEAD

pars <- list(
  t1 = -60,
  dt = 0,
  r1 = 0.04,
  r2 = 0
=======

pars <- c3pars(
  t1 = -60,
  dt2 = 0,
  k1 = 0.04,
  k2 = 0
>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba
)
tchr <- ChromCom3(pars, time = seq(from=-140, to=90, by=0.1))
tchr <- generateCells(tchr, nsim=10000)
x <- seq(from=pars$t1, to=100, by=1)
<<<<<<< HEAD
y <- exp(-pars$r1 * (x - pars$t1))
=======
y <- exp(-pars$k1 * (x - pars$t1))
>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba
g <- plotTimelines(tchr)
g + geom_line(data=data.frame(x=x, y=y), aes(x,y), colour="yellow")


<<<<<<< HEAD

=======

>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba




Parameter grid



parameterGrid <- function(pars, par1, range1, par2, range2) {
  melts <- NULL
  for(p1 in range1) {
    pars[[par1]] <- p1
    for(p2 in range2) {
      pars[[par2]] <- p2
      chr <- ChromCom3(pars)
      chr <- generateCells(chr)
      label1 <- sprintf("%s = %.3g", par1, p1)
      label2 <- sprintf("%s = %.3g", par2, p2)
      m <- meltTimelines(chr, label1=label1, label2=label2)
      melts <- rbind(melts, m)
    }
  }
  timelinePanel(melts)
}






pars <- c3pars(
  t1 = -30,
  dt2 = 0,
  k1 = 0.05,
  k2 = 0.01
)
range1 <- c(0.01, 0.05, 0.08, 0.12)
range2 <- c(0.01, 0.05, 0.08, 0.12)
parameterGrid(pars, "k1", range1, "k2", range2)


<<<<<<< HEAD

=======

>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba



I think we need to introduce a delay between pink and red.



pars <- c3pars(
  t1 = -30,
  dt2 = 10,
  k1 = 0.05,
  k2 = 0.01
)
range1 <- c(0.01, 0.03, 0.05, 0.10)
range2 <- c(0.01, 0.05, 0.08, 0.12)
parameterGrid(pars, "k1", range1, "k2", range2)


<<<<<<< HEAD

=======

>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba



What about range of delays?



pars <- c3pars(
  t1 = -30,
  dt2 = 0,
  k1 = 0.03,
  k2 = 0.01
)
range1 <- c(0, 10, 20, 30)
range2 <- c(0.01, 0.05, 0.08, 0.12)
parameterGrid(pars, "dt", range1, "k2", range2)


<<<<<<< HEAD

=======

>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba



Experimental data.



echr <- experimentalData(dataFile$scramble)
plotTimelines(echr, smooth=TRUE, k=15)


<<<<<<< HEAD

=======

>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba

13/07/2017

Prepared shiny app for deployment.

24/07/2017

Tomo suggested that the model parameters can be derived directly from data. We can write the following equations

\(\dot{B} = -k_1 B\)

\(\dot{P} = k_1 B - k_2 R\)

\(\dot{R} = k_2 P\)

\(B + P + R = 1\)

<<<<<<< HEAD

However, I’m not sure if this is exacly the model I use in the simulation. In the simulation, we assumed that P->R transition can happen only after the BB->P transition.

=======

However, I’m not sure if this is exactly the model I use in the simulation. In the simulation, we assumed that P->R transition can happen only after the B->P transition.

>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba

Because proportions and their derivatives are observed, we can find the rates:

\(k_1 = - \frac{\dot{B}}{B}\)

\(k_2 = \frac{\dot{R}}{P}\)

This will require a lot of smoothing, otherwise the derivatives will be all over the place.

First, I see if I can recover rates from the generated data.

<<<<<<< HEAD =======

timeDeriv <- function(chr, k=20) {
  ts <- list()
  for(col in chr$colours) {
    x <- chr$cnt[[col]] / chr$cnt$total
    x <- caTools::runmean(x, k)
    x <- ts(x, start=chr$timepars$start, deltat=chr$timepars$step)
    xdot <- diff(x)
    ts[[col]] <- x
    ts[[paste0(col, ".diff")]] <- xdot
  }
  ts
}

plotRates <- function(chr, k=20) {
  ts <- timeDeriv(chr, k=k)
  k1 <- -ts$B.diff / ts$B
  k2 <- ts$R.diff / ts$P
  k1[which(is.nan(k1))] <- NA
  k2[which(is.nan(k2))] <- NA
  df <- data.frame(t=time(k1), k1=k1, k2=k2)
  m <- melt(df, measure.vars = c("k1", "k2"))
  m$value <- as.numeric(m$value)
  m$value[which(m$value > 0.5)] <- 0.5
  m$value[which(m$value < -0.5)] <- -0.5
  ggplot(m, aes(x=t, y=value)) + geom_line(aes(colour=variable)) 
}

Start with a model with 10,000 cells and 0.1 min time step. In this model \(k_1 = 0.04\) and \(k_2 = 0\). Top panels show original data, lower panels - smoothed data.

grid.arrange(plotTimelines(tchr), plotRates(tchr, k=1), ncol=2)

grid.arrange(plotTimelines(tchr, smooth=TRUE, k=200), plotRates(tchr, k=200), ncol=2)

What happens if I use only 100 cells and 1 min time step in the model? This time \(k_2 = 0.04\).

chr <- generateCells(chr, nsim=100)
grid.arrange(plotTimelines(chr, k=1), plotRates(chr, k=1), ncol=2)

grid.arrange(plotTimelines(chr, smooth=TRUE, k=20), plotRates(chr, k=20), ncol=2)

And now, for real data. Here are scramble data smoothed with running mean with window size of 10 time points.

grid.arrange(plotTimelines(echr, smooth=TRUE, k=10), plotRates(echr, k=10), ncol=2)

Now, smoothing with 80 points.

grid.arrange(plotTimelines(echr, smooth=TRUE, k=80), plotRates(echr, k=80), ncol=2)

25/07/2017

Had a chat with Tomo and John yesterday. They wanted to add a score showing how far the model is from the experimental data. I did it by calculating \(\chi^2\) summed over all three curves. It is now added to the figure in the shiny app.

The next thing is to add transition P->B, with rate \(k_3\). It can happen only after time \(\Delta t_3\).

I added a simulation method (step-by-step). Now testing it.

pars <- c3pars(
  t1 = -30,
  k1 = 0.04,
  k2 = 0.04,
  k3 = 0.04,
  dt2 = 10,
  dt3 = 50
)
chr <- ChromCom3(pars)
chr <- generateCells(chr, nsim=1000, method="simulation")
plotTimelines(chr)

26/07/2017

What about fitting model to the data? “nlm” doesn’t work because model is stochastic. Method “SANN” from optim (that is stochastic simulated annealing) took very long time and went into a false minimum (negative \(k_2\)), as parameters cannot be constrained.

Finally I used “L-BFGS-B” method which allows box constraints. Alas, this method is not very good at finding the real minimum, often ending up in a local one. I’m guessing this is because of my stochastic model. Therefore, I need to run it several times and search for the best minimum. Using “mclapply” to speed it up.

I also changed the error score. \(\chi^2\) doesn’t work very well when expected counts are zero. And I cannot simply exclude them, because zero is an important number. Instead I calculate a simple RMS.

pars <- c3pars()
chr <- cacheData("fit_scramble_n3000_3par", fitChr, echr, pars, npar=3, nsim=3000, ntry=16, binDir=binDir)
plotTimelines(chr, expdata = echr, withpars=TRUE)

Now, the same, but with extra parameter, \(\Delta t_2\).

pars <- c3pars(dt2=5)
chr <- cacheData("fit_scramble_n3000_4par", fitChr, echr, pars, npar=4, nsim=3000, ntry=16, binDir=binDir)
plotTimelines(chr, expdata = echr, withpars=TRUE)

And now the other two data sets.

echr2 <- experimentalData(dataFile$NCAPD2)
pars <- c3pars()
chr2 <- cacheData("fit_NCAPD2_n3000_3par", fitChr, echr2, pars, npar=3, nsim=3000, ntry=16, ncores=8, binDir=binDir)
plotTimelines(chr2, expdata = echr2, withpars=TRUE)

pars <- c3pars(dt2=10)
chr2 <- cacheData("fit_NCAPD2_n3000_4par", fitChr, echr2, pars, npar=4, nsim=3000, ntry=16, ncores=8, binDir=binDir)
plotTimelines(chr2, expdata = echr2, withpars=TRUE)

echr3 <- experimentalData(dataFile$NCAPD3)
pars <- c3pars()
chr3 <- cacheData("fit_NCAPD3_n3000_3par", fitChr, echr3, pars, npar=3, nsim=3000, ntry=16, ncores=8, binDir=binDir)
plotTimelines(chr3, expdata = echr3, withpars=TRUE)

pars <- c3pars(dt2=10)
chr3 <- cacheData("fit_NCAPD3_n3000_4par", fitChr, echr3, pars, npar=4, nsim=3000, ntry=16, ncores=8, binDir=binDir)
plotTimelines(chr3, expdata = echr3, withpars=TRUE)

>>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba

timeDeriv <- function(chr, k=20) {
  ts <- list()
  for(col in chr$colours) {
    x <- chr$cnt[[col]] / chr$cnt$total
    x <- caTools::runmean(x, k)
    x <- ts(x, start=chr$timepars$start, deltat=chr$timepars$step)
    xdot <- diff(x)
    ts[[col]] <- x
    ts[[paste0(col, ".diff")]] <- xdot
  }
  ts
}

plotRates <- function(chr, k=20) {
  ts <- timeDeriv(chr, k=k)
  k1 <- -ts$BB.diff / ts$BB
  k2 <- ts$R.diff / ts$P
  k1[which(is.nan(k1))] <- NA
  k2[which(is.nan(k2))] <- NA
  df <- data.frame(t=time(k1), k1=k1, k2=k2)
  m <- melt(df, measure.vars = c("k1", "k2"))
  m$value <- as.numeric(m$value)
  m$value[which(m$value > 0.5)] <- 0.5
  m$value[which(m$value < -0.5)] <- -0.5
  ggplot(m, aes(x=t, y=value)) + geom_line(aes(colour=variable)) 
}

Start with a model wiht 10,000 cells and 0.1 min time step. In this model \(k_1 = 0.04\) and \(k_2 = 0\). Top panels show original data, lower panels - smoothed data.

grid.arrange(plotTimelines(tchr), plotRates(tchr, k=1), ncol=2)

grid.arrange(plotTimelines(tchr, smooth=TRUE, k=200), plotRates(tchr, k=200), ncol=2)

What happens if I use only 100 cells and 1 min time step in the model? This time \(k_2 = 0.04\).

chr <- generateCells(chr, nsim=100)
grid.arrange(plotTimelines(chr, k=1), plotRates(chr, k=1), ncol=2)

grid.arrange(plotTimelines(chr, smooth=TRUE, k=20), plotRates(chr, k=20), ncol=2)

And now, for real data. Here are scramble data smoothed with running mean with window size of 10 time points.

grid.arrange(plotTimelines(echr, smooth=TRUE, k=10), plotRates(echr, k=10), ncol=2)

Now, smoothing with 80 points.

grid.arrange(plotTimelines(echr, smooth=TRUE, k=80), plotRates(echr, k=80), ncol=2)

<<<<<<< HEAD ======= >>>>>>> 01b0e8d10be96d37a253e2c8c213bf6e968042ba